Kindara dataset

Detecting Seasonality in individual indicators

#load(paste0(IO$output_clue,"indiv_indicators.Rdata"), verbose = TRUE)
#colnames(indiv_indicators)

Detecting Seasonality in population indicators

#load(paste0(IO$output_kindara,"pop_indicators.Rdata"), verbose =  TRUE)
load(paste0(IO$tmp_kindara, "sk1.Rdata"), verbose = TRUE)
## Loading objects:
##   sk1
g = ggplot(sk1, aes(x = date, y = tot_sex)) + geom_line()
g

s = reshape(sk1, varying = list(4:ncol(sk1)), idvar = c("date","n_cycles","n_obs"), direction = "long")
s$time = colnames(sk1)[4:ncol(sk1)][s$time]
colnames(s) = c("date","n_cycles","n_obs","indicator","value")


g = ggplot(s[grep("sex",s$indicator),], aes(x = date, y = value/n_cycles, col = indicator))
g + geom_line() + ggtitle("Sexual frequency (#/number of cycles)")

g = ggplot(s[grep("sex",s$indicator),], aes(x = date, y = value/n_obs, col = indicator))
g + geom_line() + ggtitle("Sexual frequency (#/number of observations on that day)")

s$n_indicator  = s$value/s$n_obs
s$n_indicator  = s$value/s$n_cycles


signal = s$n_indicator[s$indicator == "tot_sex"]
time = par$time_num - 1
plot(time, signal, type = "l")

seasonal_analysis = detect_seasonal_pattern(time = time, signal = signal, title = "tot_sex Kindara", remove_weekly_patterns = FALSE)

seasonal_analysis = detect_seasonal_pattern(time = time, signal = signal, title = "tot_sex Kindara", remove_weekly_patterns = TRUE)

seasonal_analysis$time_series$date = par$date_seq - 365

save(seasonal_analysis, file = paste0(IO$tmp_kindara, "seasonal_analysis_tot_sex.Rdata"))
### TRYING WITH STL (Seasonal Decomposition of Time Series)

## weekly + yearly

# weekly period
signal_ts = ts(signal, frequency = 7) 
stl_w = stl(signal_ts, s.window = "periodic")
plot(stl_w)

# yearly period
signal_ts_no_weekly = ts(apply((stl_w$time.series[,2:3]),1,sum), frequency = length(signal)/3) 
stl_y = stl(signal_ts_no_weekly, s.window = "periodic")
plot(stl_y)

# bi-annual period
signal_ts_no_weekly = ts(apply((stl_w$time.series[,2:3]),1,sum), frequency = length(signal)/6) 
stl_ba = stl(signal_ts_no_weekly, s.window = "periodic")
plot(stl_ba)

##  "detrend" + weekly + yearly

# yearly period >> to detrend
signal_ts = ts(signal, frequency = length(signal)/3) 
stl_d = stl(signal_ts, s.window = "periodic")
plot(stl_d)

# weekly period
signal_detrended = ts(apply((stl_d$time.series[,2:3]),1,sum), frequency = 7) 
stl_w = stl(signal_detrended, s.window = "periodic")
plot(stl_w)

# yearly period
signal_ts_no_weekly = ts(apply((stl_w$time.series[,2:3]),1,sum), frequency = length(signal)/3) 
stl_y = stl(signal_ts_no_weekly, s.window = "periodic")
plot(stl_y)

### TRYING WITH Decompose

signal_ts = ts(signal, frequency =  length(signal)/3)
decompose_y = decompose(signal_ts)
plot(decompose_y)

# yearly
signal_ts = ts(signal, frequency =  length(signal)/3)
plot(signal_ts)

y = findfrequency(signal_ts)
y
## [1] 1
# weekly
signal_ts = ts(signal, frequency =  7)
plot(signal_ts)

y = findfrequency(signal_ts)
y
## [1] 1
# gaussian noise (negative example check)
signal_ts = ts(rnorm(length(signal)), frequency =  7)
plot(signal_ts)

y = findfrequency(signal_ts)
y
## [1] 1
# cosine of period 7 (positive example check)
signal_ts = ts(cos((1:length(signal))*2*pi/7), frequency =  7)
plot(signal_ts)

y = findfrequency(signal_ts)
y
## [1] 7
# cosine of period 7 (positive example check with noise)
signal_ts = ts(cos((1:length(signal))*2*pi/7)+rnorm(length(signal))/2, frequency =  7)
plot(signal_ts)

y = findfrequency(signal_ts)
y
## [1] 7
# trying to remove the trend first to check that it's not just the first component that predominates
seasonal_analysis = detect_seasonal_pattern(time = time, signal = signal, title = "tot_sex Kindara", remove_weekly_patterns = TRUE)

detrended_signal = signal -  seasonal_analysis$time_series$trend

# yearly
detrended_signal_ts = ts(detrended_signal, frequency =  length(signal)/3)
plot(detrended_signal_ts)

y = findfrequency(detrended_signal_ts)
y
## [1] 1
# weekly
detrended_signal_ts = ts(detrended_signal, frequency =  7)
plot(detrended_signal_ts)

y = findfrequency(detrended_signal_ts)
y
## [1] 1
# does not work either
fft_sex = fft(signal)
plot(abs(fft_sex), type = "l")

fft_sex = fft(detrended_signal)
plot(abs(fft_sex), type = "l")

period_in_days = length(signal)/(1:(length(signal)/2)) 
plot(period_in_days,abs(fft_sex)[2:(length(signal)/2+1)], type = "l", log = "x")